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ABSTRACT 

Estimates of the dynamical surface mass density at the solar Galactocentric distance are commonly 
derived assuming that the disk is in vertical equilibrium with the Galactic potential. This assumption 
has recently been called into question, based on the claim that the ratio between the kinetic and the 
gravitational energy in such solutions is a factor of 3 larger than required if Virial equilibrium is to hold. 
Here we show that this ratio between energies was overerestimated and that the disk solutions are likely 
to be in Virial equilibrium after all. We additionally demonstrate, using one-dimensional numerical 
simulations, that the disks are indeed in equilibrium. Hence, given the uncertainties, we find no reason 
to cast doubt on the steady-state solutions which are traditionally used to measure the matter density 
of the disk. 

Subject headings: Galaxy: disk - Galaxy: kinematics and dynamics - Galaxy: structure 



1. INTRODUCTION 

Studies on the amount of unseen mass associated with 
the Galactic disk, have been carried out since Oort (1932) 
and up to the present day (e.g., Kuijken & Gilmore 1991; 
Bahcah et al. 1992; Creze et al. 1998; Siebert et al. 2003; 
Holmberg & Flynn 2004; Kalberla et al. 2007; Moni Bidin 
et al. 2010). All recent studies indicate that there is no 
dynamically significant amount of disk dark matter at the 
solar position, although a small fraction (<10%) cannot be 
ruled out. In order to derive the dynamical mass surface 
density, the usual way to proceed is to consider a tracer of 
the potential, such as K giant stars or Hi gas. In princi- 
ple, ideal tracers are old enough (age > 1 Gyr) that they 
are expected to be well mixed into the Galactic potential, 
so that they satisfy the Virial condition. In a recent pa- 
per, Garrido Pestaiia & Eckhardt (2010) (hereafter GPE) 
compare the kinetic energy T of stellar and gaseous com- 
ponents in the local Galactic disk, to their gravitational 
energy VF, and conclude that the vertical configuration is 
not in a state of equilibrium because |2T/IF| ~ 3 and, 
therefore, the classical Poisson-Boltzmann approach used 
to interpret the data is an inappropriate and futile exer- 
cise. If the energy in vertical random motions is a factor 
of 3 larger than in Virial equilibrium, the local portion of 
the disk should be in a phase of violent relaxation. Un- 
der such conditions, the thickness of the disk would inflate 
by a factor of 5/3 on average (see section [5]) after just a 
few vertical crossing times (~ 10* years). It is unlikely 
that we would be observing the disk at a time just prior 
to such a violent event, and if we were, the classical view 
that the disk thickening is due to slow dynamical heat- 
ing due to encounters with massive molecular clouds, spi- 
ral arms and gravitational perturbations caused by minor 
mergers, should change dramatically. In this paper, we 

^ Note that W > in this case. GPE uses a different sign convention. 



show, however, that there is no reason for abandoning the 
assumption that the vertical configuration is in a steady 
state. 

2. STEADY-STATE CONFIGURATION 

Consider a rotationally symmetric disk stratified in the 
z-direction, where z = is the disk plane of symmetry, and 
composed of N components in equilibrium. Each compo- 
nent, either if it is collisionless (stars and dark matter) 
or collisional (gas), satisfies a Boltzmann equation. The 
vertical Jeans equations can be combined to obtain: 
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dz 
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where iJi{z) is the vertical velocity dispersion of the com- 
ponent i, ^{z) is the gravitational potential and p{z) = 
'^Pi{z). For gaseous components, Oi is the effective ve- 
locity dispersion, which includes the pressure support by 
turbulent motions, magnetic fields and cosmic rays. In a 
real galaxy, gradients in the radial and azimuthal direc- 
tions generate nondiagonal terms of the velocity disper- 
sion tensor, i.e. cr/jz and ugz-, in the vertical Jeans equa- 
tion. At z < 1 kpc, the fractional corrections due to the 
cross terms are less than 1% and are usually ignored ex- 
cept when studying stars at z > 2 kpc (e.g., Bahcall 1984; 
Oiling 1995; Becquaert & Combes 1997; Kalberla 2003; 
Moni Bidin et al. 2010). 

Following GPE, we start by considering a infinite plane- 
parallel layer of total column density Sqo in equilibrium 
with its own gravity. By choosing the zero of the potential 
at z = and multiplying the above equation by z and inte- 
grating over z, one obtains the Virial theorem, linking the 
kinetic and gravitational energies with surface terms. If 
the surface terms are null at z — >■ ±oo, the relationship is 
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2T = W, where W and T are the gravitational and kinetic 
energies (due to vertical motions) , respectively, per surface 
unit in the dislfl- GPE show that the Virial theorem can 
be rewritten as R = Q where 



R 



and 



ap{z)p(z + a) dz da, 



(2) 
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where p{z) = p{z) /T,oo- Here (ct^) is the mass-weighted 
mean-square velocity dispersion, i.e. (c^) = 2T/Eoo- 

For the particular case of having a single self-gravitating 
isothermal component {N = 1), it holds that Q is exactly 
1/2, whereas it is slightly different from 1/2 for a Gaussian 
density profile in equilibrium {Q = 0.45). GPE argue that 
Q ~ 1/2 and evaluate the left-hand-side of Eq. 1^ for the 
Galactic disk at the solar position and find that it has the 
value 1.6. Since this is a factor of ^ 3 larger than required 
if the Virial equilibrium is to hold, they conclude that the 
disk at the solar position is not in a steady-state vertical 
configuratior0. 

In the following, we will check the robustness of the as- 
sumption that Q ~ 1/2. To do so, we start by integrating 
Equation ([Ij from the midplane to a certain distance z 
above the midplane and obtain: 

N 

1=1 

For a system such that the volume density of kinetic en- 
ergy decays to zero at large z, Eq. (U) becomes 
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and tell us that the total pressure at the midplane must be 
in balance with the weight of all the material above z — Q. 

In the particular case of a stratified plane-parallel layer, 
all the variables depend only on z. Integration of the Pois- 
son equation allows us to find the vertical acceleration. If 
the layer is confined by its own gravity, the vertical accel- 
eration for z > reads 
d<i> 

27rGS(z), (6) 



dz 

where S(z) is defined for z > as 
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Combining Eqs. ([S]) and and using the relation 
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we find 
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Note that the above equation is exact for a plane-parallel 
layer in equilibrium. Therefore, the total pressure at z = 
(not p(0) (cr2) /(ttGS^) as stated in GPE) divided by 
ttGS^ must be 1/2 in a steady-state plane-parallel layer. 
GPE analysis applied to a plane-parallel layer is strictly 
correct only if p(0) (cr^) = Pj(0)o'i^(0). We will see later 
that this equality does not hold in a multicomponent disk. 

In a real galaxy, not only the disk but also other com- 
ponents (such as the dark matter halo) contribute to the 
vertical confinement of stars in the thin and thick disks. 
From the Poisson equation, the gravitational vertical ac- 
celeration is given by 



^ = 2TiG[^d{z) + ^h{z)] 
oz 



n{n + G)dz\ (11) 



where Sd(-z) and Sh(z) are the column density of the 
disk and of a spherical dark halo, respectively, = 
R-^d^/dR and G = R{dVL/dR) (e.g., Kuijken & Gilmore 
1989; Ferriere 1998). At i? > 5 kpc and within a few kpc 
from the disk plane, -f G ~ (e.g., Kuijken & Gilmore 
1989) and thus 



(9$ 



= 27rGS(z) 
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dz ' ' R Vi?2 + z2 ' 

where Vfi{R) is the contribution of the dark halo to the ro- 
tation curve; Vh — 140 km s^^ at the solar position (e.g., 
Flynn et al. 1996; Klypin et al. 2002). In the particular 
case of a massless disk embbeded in the spherical poten- 
tial of the isothermal dark halo, we recover the classical 
result that the vertical epicyclic frequency coincides with 
the orbital frequency at the midplane, d'^^/dz'^\z=() — 
Substituting Eq. (fT2|) into Eq. © and for z <^ R, we ob- 
tain 



EpKo)^?(o) 
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which can be expressed as 

ttGJ^I, 2- 

^ If true, the disk will reach a new (final) equilibrium configuration, with kinetic energy Tf and potential energy Wj, such that E = (3/2)Wj, 
where E is the kinetic plus gravitational energies. In the initial nonequilibrium configuration, E = Ti + Wi = (5/2)Wi. Equating the energy, 
we find that Wf = {h/'i)Wi. Since the scaleheight scales roughly as W , the disk will become thicker by a factor of 5/3. 



where the index i refers to the disk components, hi is the 
vertical scaleheight of population i and = Sd(z = oo). 

Evaluation of X] Pii^)'^} (not including the stellar halo) 
in the updated disk model of Holmberg & Flynn (2004) re- 
ported in Table 2 of Flynn et al. (2006) gives 20.2Mopc-3 
km^ s^^. Taking the total surface density in disk compo- 
nents of Soo = 48.7M0pc~^, the left-hand-side of Equa- 
tion (fT3|) is 0.61. We then compute the right-hand-side of 
Eq. ([13]) at i? = i?0 = 8 kpc and find that it is ~ 0.56, 
which is only 9% smaller than the left-hand-side. This 
simple estimate indicates that the disk model of Flynn et 
al. (2006) is very close to satisfying the Virial theorem. In 
order to assess if this difference of ten percent is significant, 
we have calculated Pi{z) in the model of Flynn et al. (2006) 
and computed the right-hand-side of Eq. (fT5|) exactly. We 
confirmed that the disk model fulfills the Virial condition. 

Why did GPE obtain that the disk models by Bahcall 
et al. (1992) and Flynn et al. (2006) do not comply with 
the Virial theorem? It was a consequence of the combina- 
tion of three factors: (1) they adopted Q = 1/2, which is 
in error by a factor of 2 (see below), (2) the stellar halo 
was treated as a disk component, which introduces an ad- 
ditional artificial factor of 1.25 and (3) the contribution of 
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the dark halo in the vertical confinement of disk stars was 
not taken into account. As a consequence, the ratio be- 
tween kinetic and gravitational energies was overestimated 
by a factor of 3. 

For a stratified plane-parallel self-gravitating layer, it is 
true that i? = Q in Virial equilibrium. However, Q is not 
necessarily 1/2 in a generic disk. In fact, the form of the 
left-hand-side of the equilibrium condition used in GPE 
(Eq. [2]) can be preserved by introducing a "correction" 
factor e on the right-hand-side: 

P(0)(^^) _ e 
ttGE^ 2' 

with 



(14) 



p(0)(a2) 



(15) 



Eft (OK'' 

The value of e depends on the vertical structure of the disk. 
For a single Gaussian (non-isothermal) density profile in 
z, e = (cr^) /ct2(0) = 2\/2/7r = 0.9, and hence Q = 0.45, 
as reported in GPE. However, for the 14th-components 
of the updated disk model of Holmberg & Flynn (2004), 
the e-value in the Virial theorem criterion for a state of 
equilibrium is 2 and thus Q ~ 1. To be more precise, we 
adopted the values of pi (0) and ct; for the 14th-components 
of Flynn et al. (2006) and solved the Poisson-Boltzmann 
equation for such a disk with no dark halo, in order to com- 
pare with GPE analysis on a common ground. We obtain 
p(0) = O.O91A/r0pc-3, = 55Mopc-2, {a^f'"^ = 22.25 

km s^^ and R = Q = 1.1. Therefore, GPE analysis applied 
to the infinite plane-parallel multicomponent layer overes- 
timates the T/W ratio by a factor of 2.2 because 



= 2.2. 



3. SIMULATIONS 

In the last section we have shown that the Virial theorem 
is satisfied for the Galactic disk models under considera- 
tion. However, the Virial theorem is a relation between 
global quantities and is not a guarantee that the system 
is in a steady-state equilibrium. In order to dispel any 
doubt, we have tested the temporal stability of the type of 
disks considered here using N-body simulations, in which 
stars are initially set up in a one-dimensional density and 
velocity distribution, based on solutions to the Poisson- 
Boltzmann equation, after which we calculate their verti- 
cal motions under the gravitational force generated by the 
overall density distribution of the disk. 

At each time step in the simulation, the downward accel- 
eration at height z above the disk, is found by integrating 
numerically the density profile of all the stars (Eq. ^ to 
absolute height \z\. A leap-frog integrator is then used to 
compute new positions and velocities for the particles, and 
the process repeated. We keep track of the velocity disper- 
sion, scale height and other quantities for the ensemble, to 
check if it is in vertical equilibrium or not. 

Three disk types were considered. Firstly, we looked at 
an "isothermal disk", in which the initial conditions are 
a sech^(z//i) distribution in density (with a scale height 
h) and a (constant) velocity dispersion cr(z), given by 
by cr = h^2TiGp{Q), where p{Q) is the disk central den- 
sity. The density distribution in such a disk is given by 



p{z) = p(0)sech2(z/2/i). We adopt p(0) = 0.1 Mq pc'^ 
and h = 250 pc, implying a total disk surface density of 
50 M0 pc~^ and a velocity dispersion of cr = 13.0 km 
s~^. These are similar to the measured values for the local 
Galactic disk. 

Secondly, we looked at an "exponential disk", in 
which the density distribution is given by p{z) = 
p(0)exp(— |z|/ft,), and the velocity dispersion is a function 
of z, being smaller close to the midplane and rising asymp- 
totically to a constant value at large height. The velocity 
dispersion a{z) is given by 



cr(z) = 2hy/TrGp{0){l - 0.5exp(-|z|//i)). (17) 

Adopting a scaleheight of h — 250 pc and a central den- 
sity of O.IM0 pc~'^, this disk has a surface density of 50 
Mq pc~^, like the isothermal disk. The central velocity 
dispersion for this disk rises from 13.0 km s^^ at z = to 
18.4 km s~^ at large z. 

The simulations of these two disks were run with 10 
million particles, for 5 Gyr, over 1 Myr time steps. This 
time step was found to be a good compromise between the 
length of the runs and the accuracy of the orbits. Note 
that the oscillation time for stars in the tested potentials 
is typically 30 to 80 Myr. Both disks are completely self- 
gravitating, with the gravitational field provided by the 
density distribution of the particles alone. 

Over 5 Gyr, the central density, mean velocity disper- 
sion and scaleheights of the simulated disks remain within 
a few percent of their initial values. Thus, the exponential 
disk and the isothermal disk are found to be stationary in 
the simulations. Furthermore, this is what we expect, since 
GPE's analysis shows that, for both simulated disks, the 
Virial equilibrium condition is met very precisely. Thus, 
no long term secular evolution is expected for either disk, 
and this we fully confirm with the simulations. 

We then simulated the much more complex system, con- 
sisting of all components, stellar and gaseous, of the disk 
in Flynn et al. (2006). This disk consists of 14 isothermal 
components, with a range of central densities and (isother- 
mal) velocity distributions. The stellar halo, which is the 
15th component of the Flynn et al. (2006) disk model, is 
ignored (including it in the solutions is found to have no 
effect on the disk stability, but we focus in the following 
on disk components only). The gaseous components are 
modelled as coUisionless particles, just like the stars, since 
the aim is to test the Poisson-Boltzmann solution, which 
implicitly makes this assumption. We also include the ef- 
fect of the halo dark matter in these simulations, not by 
simulating it with particles, but by adding it in to the 
computation of the vertical force as an extra density term 
(as a function of z). 

The simulation is run with 10 million particles for a pe- 
riod of 5 Gyr with time steps of 1 Myr. Energy conserva- 
tion in the simulation is better than 0.1 percent. We find 
that this complex disk is also stationary. The density and 
velocity distributions of the particles after 5 Gyr are very 
similar to their initial distributions, and so no large scale 
secular evolution takes place. If this rather complex disk 
were substantially out of equilibrium initially, we would 
expect rapid evolution in these quantities. Figure [1] shows 
the results. Panel (a) shows the evolution of kinetic en- 
ergy (dotted line), potential energy (dashed line) and to- 
tal energy (solid line) over the 5 Gyr of the simulation. 
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Panel (b) shows the total velocity dispersion as a function 
of time. In panels (c) and (d), we show the initial (cir- 
cles) and final (crosses) density and velocity distributions 
of the particles as a function of z height. The initial and 
final distributions are essentially indistinguishable. Panel 
(e) shows the time evolution of the central density, and 
panel (f) the scale height (/i, of an exponential fit to the 
density profile, measured from the midplane to z « 750 
pc). None of the plots show any significant evolution in 
these quantities with time. We conclude from the simula- 
tions that the Flynn et al. (2006) model disk is indeed in 
equilibrium, as we concluded from theoretical analysis in 
the previous section, and contrary to the claims in GPE. 

In order to illustrate the role of the dark halo in con- 
fining the disk, we rerun the simulation but the dark halo 
contribution to the force was turned off so that the disk is 
not in equilibrium initially. Figure [5] shows that the disk 
relaxes dynamically (particularly those stars in the thick 
disk) to a stationary state in 0.3 Gyr. Panels (a) and (b) 
show the typical behaviour of T, ly and {<J^) in a relax- 
ing process. As can be seen from panels (c), (d) and (e), 
the disk becomes thicker and colder. Panel (f) shows the 
temporal evolution of R and Q. Since all the gravity is 
provided by the stars and the gas only, Virial theorem for 
this system dictates that R = Q. We see that R = Q — I 
in the steady-state. 



4. DISCUSSION AND CONCLUSIONS 

We have examined the claim by GPE that solutions for 
the density distribution of matter in the local Galactic 
disk are seriously in error because such solutions are not 
in Virial equilibrium. GPE claim that the ratio between 
the kinetic and the gravitational energy in such disks is a 
factor of 3 larger than required for Virial equilibrium to 
hold. We have shown that this ratio between energies was 
overerestimated by GPE because (1) they assumed that 
Q — 1/2, (2) the stellar halo was treated as a disk com- 
ponent and (3) the vertical confinement provided by the 
dark halo was ignored. 

Taking into account these corrections, the local disk 
model of Flynn et al. (2006) is likely to be in Virial equlib- 
rium after all. Our technique could easily be applied to any 
such model (e.g., Bahcall et al. 1992; Just & Jahreiss 2010) 
and yield a similar result. Simulations of the oscillation of 
stars in model disks show that the initial distribution of the 
stars in density and velocity space are stationary, backing 
up our theoretical analysis. 

We would like to thank J. L. Garrido and the referee for 
very enlightening comments. CF is very grateful to the 
University of Sydney, for hosting a visit where this work 
was carried out. This work has been partly supported by 
CONACyT project CB-60526F. 
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Fig. 1. — Simulation of the secular evolution of the Flynn et al. (2006) model disk. Panel (a) shows the total energy (solid line), the kinetic 
energy (dotted) and the potential energy (dashed) as functions of time, from to 5 Gyr. Panel (b) shows the total velocity dispersion of all 
the particles as a function of time. The central panels, (c) and (d), show the initial (circles) and final (crosses) density and velocity dispersion 
distributions of the particles as functions of height above the disk z; the distributions are close to identical, indicating that the disk solutions 
are stationary and that there is no long term evolution in its structure. Panels (c) and (f) show the evolution of central density and scale 
height as functions of time. All indicators demonstrate that the disk is in a stationary state in the long term. 
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Fig. 2. — Simulation of the secular evolution of tiie Flynn et al. (2006) model disk, as in Figurc[T] but with the special case that the dark 
matter halo has been removed from the simulation. All panels are as for Figure [T] except panel (f), which shows the quantities R (dashed 
line) and Q (solid line). Note that the time axis shows to 1 Gyr, to highlight the rapid response of the disk. The disk is initially out of 
Virial equilibrium, but adjusts itself to equilibrium within a few dynamical times, arriving at the state R = Q. Since the gravity of the dark 
halo is lacking, which would help to confine the disk, the response is for the potential energy to increase which kinetic energy decreases, i.e. 
the disk thickens (panel c) and becomes kinematically cooler (panel d). 



